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Abstract 

The steady flow of three viscoelastic fluids (Oldroyd-B, FENE-P, and Owens model 
for blood) in a two-dimensional chaimel, partly bound by a deformable, finite 
thickness neo-Hookean solid, is computed. The limiting Weissenberg number be¬ 
yond which computations fail to converge is found to increase with increasing 
dimensionless solid elasticity parameter E, following the trend Owens > FENE- 
P > Oldroyd-B. The highly shear thiiming nature of Owens model leads to the 
elastic solid always collapsing into the charmel, for the wide range of values of E 
considered here. In the case of the FENE-P and Oldroyd-B models, however, the 
fluid-solid interface can be either within the charmel, or bulge outwards, depending 
on the value of E. This behaviour differs considerably from predictions of earlier 
models that treat the deformable solid as a zero-thickness membrane, in which case 
the membrane always lies within the charmel. The capacity of the solid wall to 
support both pressure and shear stress, in contrast to the zero-thickness membrane 
that only responds to pressure, is responsible for the observed difference. Compar¬ 
ison of the stress and velocity fields in the channel for the three viscoelastic fluids, 
with the predictions for a Newtonian fluid, reveals that shear thinning rather than 
elasticity is the key source of the observed differences in behaviour. 
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1. Introduction 

Numerical simulation of blood transportation through the human cardiovas¬ 
cular system is an intense area of research M. Blood, which is rheologically 
complex, interacts with blood vessels walls both chemically and mechanically to 
give rise to an intricate fluid-structure interaction. From a fluid mechanics point 
of view, blood flow in large diameter blood vessels is commonly referred to as the 
macrocirculation, while flow in small vessels, such as arterioles, venules, and capil¬ 
laries is referred to as the microcirculation. The Navier-Stokes equations are a good 
model for blood flow in the medium to large arteries, since Reynolds numbers are 
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high, tube diameters are large, and blood can be considered to be an incompressible 
viscous Newtonian fluid. Most studies in the literature on fluid-structure interac¬ 
tion in the context of blood flow have so far been focused on the macrocirculation. 
In small vessels, however, where the shear rate is small, blood behaves as a non- 
Newtonian fluid due to its particulate nature. This necessitates its shear thinning, 
viscoelastic and thixotropic nature to be taken into account 10-0]. Another impor¬ 
tant aspect of the microcirculation is that the vessel wall thickness to diameter ratio 
is very high |@]. To the best of our knowledge, so far there have been no studies of 
fluid-structure interaction associated with the flow of viscoelastic fluids in vessels 
with finite thickness walls. The aim of this work is to examine the flow of a variety 
of viscoelastic fluid models interacting with a finite-thickness elastic vessel wall. 

The study of fluid flow in collapsible charmels and tubes has been motivated by 
the complex and nonlinear dynamics revealed by laboratory experiments 
The earliest and simplest theoretical models of collapsible-tube flow were lumped- 
parameter 111] and one dimensional models Em, followed by two-dimensional 
models where part of one wall is replaced by a tensioned membrane in a two- 
dimensional rigid parallel sided charmel. The membrane model assumes that 
the bending stiffness and extensibility of the wall in the flow direction can be 
ignored, and that the movement of the elastic wall is only in the direction normal 
to the wall. More recently, this basic model has been improved by using a plane 
strained elastic beam model for the collapsible wall with a Bemoulli-Euler beam, a 
Timoshenko beam and a 2D solid model 114-1^. Wall stiffness was found to play 
a major role in attaining a steady state solution, and for very small wall stiffness, 
the results of the beam model compared favourably with those of the membrane 
model. Work is ongoing on extending these models to describe 3D compliant 
tubes |Il6-21|. In all these studies, however, the fluid has always been treated 
as Newtonian. In this work, we examine the flow in a 2D collapsible charmel of 
three different viscoelastic fluid models, namely, the Oldroyd-B, the FENE-P and 
the Owens model for blood |0|, by considering the deformable wall to be a finite¬ 
thickness incompressible neo-EIookean solid. The charmel dimensions are chosen 
to be compatible with the microcirculation. 

The steady flow of Oldroyd-B, FENE-P and the Owens model fluids in a two- 
dimensional collapsible charmel has been studied recently by Chakraborty et al. 
1220. In contrast to the present work, however, the coll aps ible wall was modelled 
as a zero thickness membrane under constant tension 12311 . It was shown that the 
predictions for the different viscoelastic fluids differ significantly from each other, 
with the key factor being the extent of shear thinning predicted by the individual 
models. In particular, it was shown that viscoelastic fluids behave identically to 
Newtonian fluids, provided that the viscosity of the two fluids at the location of 
the maximum shear rate in the charmel is the same. Subsequently, Chakraborty 
and Prakash ll23l examined the influence of the degree of shear thinning of the 
viscoelastic fluid by systematically varying the finite extensibility parameter bu, 
in the FENE-P model, which controls the extent of shear thinning experienced by 
the fluid, and is consequently a convenient parameter for examining the influence 
of shear thirming. They found that the pressure drop, the molecular conformation 
tensor fields, and the stresses in the flow domain are significantly affected by the 
extent of shear thinning of the FENE-P fluid. Importantly, in both these stud¬ 
ies, it was found that the significant differences that arise amongst the different 
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Figure 1: Geometry of the domain. = 7W, L = 5W, and Lj = 7W. The wall thickness t is varied 
between 0.1 W and 0.4W. 


viscoelastic fluids in the predicted value of the tangential shear stress on the mem¬ 
brane surface, has no influence on the shape of the deformable membrane, because 
of the boundary condition adopted in the work. Essentially it was assumed that 
the shape of the membrane is governed only by the normal stresses acting on it. 

In order to use a more realistic model for the collapsible wall, Chakraborty et al. 
replaced the zero thickness membrane with a finite thickness neo-Hookean 
solid, which can account for the effect of shear stress on the shape of the wall. 
The fluid was, however, assumed to be Newtonian. The model formulation fol¬ 
lowed the seminal work of Carvalho and Scriven 11261] . who examined roll cover 
deformation in roll coating flows, with the rubber roll cover modelled as either 
an incompressible neo-Hookean or a Mooney-Rivlin solid. Computational predic¬ 
tions of the deformation of the collapsible wall and pressure drop in the channel 
were found to be in good agreement with experimental measurements carried out 
in a polydimethylsiloxane microfluidic device, composed of a single microchaimel 
with a thin flexible layer along one side of the chaimel. 

The present work extends the model of Chakraborty et al. ll2^ by replacing the 
Newtonian fluid with Oldroyd-B, FENE-P and the Owens model viscoelastic fluids. 
The constitutive equations of all the viscoelastic fluids are written in conformation 
tensor form ll^ . and the fully coupled, steady state, fluid and solid equations 
are solved using the DEVSS-TG/SUPG (discrete elastic viscous stress split-traceless 
gradient) finite element method. Such a general formulation enables us to study 
the complex fluid-structure interaction that arises both from the capacity of the 
deformable wall to support shear and normal stresses, and from the shear thinning 
and viscoelastic nature of the fluid. 

The plan of the paper is as follows. The problem formulation, with details of the 
governing equations for the viscoelastic fluids and incompressible neo-Hookean 
solid, the boundary conditions and the relevant dimensionless variables are pre¬ 
sented in Section |2l The results of viscoelastic and Newtonian fluid computations 
are compared in Section |3]. In particular, the dependence of the shape of the fluid- 
solid interface, and of the pressure, stress, conformation tensor and velocity fields 
on the different parameters, is examined. Finally, concluding remarks are drawn 
in Section m 
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2. Problem formulation 

The geometry of the flow is that of a 2D channel, with one of the walls con¬ 
taining an elastic segment as illustrated in Fig. [U In units of chaimel width W, the 
dimensions of the chaimel are = 7W, L = 5W, and = 7W. In most of the 
simulations, the solid wall has a thickness t = 0.41V. 


2.1. Governing Equations 

We have nondimensionalized the various physical quantities, by scaling lengths 
and displacements with W, velocities with GW/qo and pressure and stresses with 
G. Here G is the shear modulus of the solid, 7]o(= r/s + '?p,o) is the zero shear rate 
solution viscosity, ijs is the solvent viscosity and 7]p,o is the contribution of the micro- 
structural elements to the zero shear rate viscosity. (For a Newtonian fluid, 7]o is 
just the constant Newtonian viscosity). Non-dimensionalization of the governing 
equations and boundary conditions yields the following dimensionless numbers: 


Re 


pWUo 

no 


2l; Wi 

m 


AqUo 

W 


no^o 

GW 


p=l- 

G 


( 1 ) 


where Re is the Reynolds number, jS is the viscosity ratio, Wz is the inlet Weissenberg 
number, F is the dimensionless solid elasticity parameter, P is the dimensionless 
pressure, p is the density of the liquid, Uq is the average inlet velocity, Aq is the 
constant characteristic relaxation timy of the microstructure. It is also convenient to 
define a local Weissenberg number Wz = Aq y, which measures the non-dimensional 
shear rate at any location in the flow. 

Upon introduction of these dimensionless variables, we can recast governing 
equations in the following dimensionless form: 


V • v = 0 (Mass balance) (2) 

Re 

— V Vv = V • T (Momentum balance) (3) 

T F 

V VM- Vv^ - M-M - Vv = {/(trM)M-I } (Conformation tensor ) (4) 

T = -PI -I- T® -I- rP (Cauchy stress tensor) (5) 

T® = jS(Vv -I- Vv^) (Viscous stress tensor) (6) 

F 

tP = (1-^) —{/(trM)M-I} (Elastic stress tensor) (7) 

Vx -8 = 0 (Equation of motion for solid) (8) 

S = F“^ • cr (First Piola-Kirchhoff stress tensor ) (9) 

cr = -Til -I- B (Cauchy stress tensor for a neo-Hookean material) (10) 

In these equations, v is the velocity, V denotes the gradient, I is the identity tensor, 
71 is the pressure, and B is the left Cauchy-Green tensor, expressed as B = F • F^. 
The deformation gradient tensor F relates the undeformed state [X = (X, Y, Z)] to 
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the deformed state [x = {x, y, z)] and is expressed as: 


F = 


dx 

dX 


( 11 ) 


The form of /(tr M) is model specific. For the Oldroyd-B model, /(tr M) = 1, 
while for the FENE-P model 


(-M- — 

where, Bm is the finite extensibility parameter, defined as the ratio of maximum 
length squared of the micro-structural element to its average length squared at 
equilibrium. 

Eor the Owens model, /(tr M) = 1 also holds, however, the constant relaxation 
time Ao in Eq. (HJ is replaced by a function A, which represents the relaxation time 
of the elastic stress due to blood cell aggregates. Note that this replacement is not 
carried out in Eq. dZl), where the relaxation time remains constant and equal to Aq. 
The function A depends on the average size of the blood cell aggregates, n, which 
is controlled by the competition of spontaneous aggregation and flow-induced 
disaggregation. lolov et al. iIa^ have recently developed a finite element method 
for solving the Owens model in its complete generality. Since our focus here is 
on developing a fluid-structure interaction model that accounts for a viscoelastic 
fluid model and a finite thickness elastic wall, we assume for simplicity that the 
dynamics of n are fast with respect to other changes of the flow, i.e., n = nst(y), 
which is its equilibrium value based on the local shear rate / = V2 D : D, where D = 
i(Vv -I- Vv^) is the rate of strain tensor. This choice preserves the viscoelastic and 
shear thiiming character of blood but does not capture its thixotropic behaviour |@]. 
This simplification makes it unnecessary to solve an additional equation for the 
variation of n in the flow domain. Under this assumption, the relaxation time A is 


A = 


A 


H 




ripif) 


(13) 


where. Ah is the relaxation time of individual blood cell aggregates, is the 
infinite shear-rate viscosity, and r]p(y) is their contribution to blood viscosity given 
by the Cross model. 


/?p(y) = ripfi 


1 + 0iy”'\ 
1 + 02 /'"/ 


(14) 


where m is a power law index, and the ratio of parameters 0i and 02 satisfies the 
expression, 0i/02 = ?]p,oo/Tip,o IZaL More details on the current implementation of 
Owens model are given in Ref. 1301]. The values of all model parameters used here 
are reported in section |3^ 

In free boundary flow problems, one of the major difficulties lies in the fact that 
the location of the free boundaries is unknown a priori and their solution is a part of 
the total solution. Here we use a boundary fitted finite element based elliptic mesh 
generation method |l3ll-l^ which involves solving the following elliptic differential 
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equation for the mapping: 


V • D • = 0 


(15) 


where ^ is a vector of positions in the computational domain and the dyadic, D, is 
a function of analogous to a diffusion coefficient, which controls the spacing of 
the coordinate lines 113311 . 


As mentioned earlier, the formulation of the fluid-structure interaction problem 
posed here follows the procedure introduced previously by Carvalho and Scriven 
j2(i\ in their examination of roll cover deformation in roll coating flows. However, 
it turns out that the weighted residual form of Eq. ((HI) used in their finite element 
formulation is incorrect. While the error does not lead to significant discrepancies 
for small deformations, it is serious for large deformations. The correct form of the 
weighted-residual equation is presented in Chakraborty et al. |j^. 

2.2. Boundary conditions and discretization 

We prescribe the following boundary conditions: 

1. No slip boundary conditions (v = 0) are applied on the rigid walls. 

2. Zero displacements are prescribed at the left side and right side of the solid. 

3. At the upstream boundary, a fully developed dimensional velocity profile is 
specified in the form, Vy = 0 and Vx = Uof{y/W) where Uq is the average 
inlet velocity. In dimensionless form this can be represented as Vx = r/(i//W) 
where, Vx is now non-dimensional. Since, for all the Wi considered here, 
the upstream velocity profiles for the Oldroyd-B and FENE-P fluids do not 
differ significantly from that of a Newtonian fluid, a Newtonian velocity 
profile is used. However, for the Owens model fluid, which is strongly shear 
thinning in nature, a different procedure is used. Since the viscosity of the 
Owens model fluid obeys power law scaling with shear rate at relatively low 
values of shear rate in simple shear flow (see section I3.1|l . we assume that 
the velocity profile at the upstream boundary is identical to the analytically 
computed velocity profile for a power-law fluid (with power law index m), 
flowing in a 2D charmel. We find in our simulations (since the entrance 
length has been assumed to be sufficiently long), that the fully developed 
velocity profile of the Owens model fluid at the inlet to the collapsible charmel 
is relatively unchanged from the analytical velocity profile imposed at the 
upstream boundary. 

4. At the downstream boundary, the fully developed flow boundary condition 
is imposed, n • Vv = 0 where n is the unit normal to the outlet. 

5. The conformation tensor equation for all the viscoelastic models is hyper¬ 

bolic in nature and therefore boundary conditions are needed only on inflow 
boundaries. At the upstream inflow, the conformation tensor does not change 
along the streamlines because the flow is fully developed Thus, 


v VM = 0 


(16) 


6. A force balance and a no-penetration condition are prescribed at the interface 
between the liquid and solid domain. 


n.T = n.cr and VsoUd = v 


fluid 


(17) 
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where n is the unit normal to the deformed solid surface. 
7. A force balance is prescribed at the top surface. 


n.cr = -PgU 


( 18 ) 


where Pf, is the dimensionless external pressure. 

8. The non-dimensional pressure of the fluid at the downstream boundary, P^, 
is set equal to zero. 


Equations (O, Q, (|1]) and (T5\i are converted into a set of algebraic equations 


by the DEVSS-TG finite element method 
interpolated velocity gradient liH 


3611 ■ which introduces the traceless 


L- 


Vv + — (V • v) I = 0 
tr 1 


(19) 


In the transport equations the rate of strain tensor D is calculated from the inter¬ 
polated velocity gradient L. 

The weighted residual form of Eqs. (O, (jl]), ((8l), ([15]) and (IT^ . yields a large 
set of coupled non-linear algebraic equations, which is solved by Newton's method 
with analytical Tacobian, frontal solver, and first order arclength continuation in 


parameters |l22l. l34l. 1371. 1381] . Note that while Eqs. (|2j), (|3|), (|4j), (ITSl) and (fl^ are 


implemented in the fluid domain, Eq. (|8|) is solved in the solid domain. 


3. Results and discussions 


A thorough validation of the finite-element code, in the context of a Newto¬ 
nian fluid, has been carried out by Chakraborty et al. 1251] . who have compared 
results obtained with the present formulation with several earlier results obtained 
in different contexts. 

It is appropriate to briefly discuss the fluid models used in the present work 
before presenting the results of our simulations, since the differences in behaviour 
amongst them is essentially due to differences in their rheology. 


3.1. Fluid models and choice of parameters 

Each of the three fluids examined here has distinct qualitative features: (i) The 
Oldroyd-B fluid is viscoelastic, but does not shear thin. Eurthermore, its uniaxial 
extensional viscosity is unbounded, (ii) The EENE-P fluid is viscoelastic, shear 
thins, and has a bounded uniaxial extensional viscosity, (iii) The Owens model 
fluid is viscoelastic and shear thins, but has an unbounded uniaxial extensional 
viscosity. Additionally, a notable feature of the Owens model, which belongs to 
the class of White-Metzner fluids, is that the dependence of viscosity rjp on shear 
rate f can be prescribed arbitrarily through the choice of parameters in the Cross 
model (see Eq. (IT4|) 1. In particular, the viscosity can be prescribed independently 
of the relaxation time. In contrast, for the EENE-P model, the dependence of rjp on 
the shear rate y is completely determined by the choice of the parameters, 7]p^o/ ^O/ 
and the finite extensibility parameter bu- Unlike in the case of the Owens model, 
no further control can be exerted on the shape of the viscosity function. 

The difference in the prediction of viscosity by the PENE-P and Owens models, 
as a function of the Weissenberg number Wi = Aoy, in steady shear flow, is displayed 
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Figure 2: The contribution of the microstructure to the total viscosity, i], for the Owens 
model and FENE-P fluids in steady shear flow as a function of Weissenberg number Wf. 
The inset shows the shear rate dependence of viscosity in the Owens model, fitted to the 
experimental results for blood reported by Chien jH], and the predictions of the EENE-P 
model for = 2 and Aq = 0.263. 


in figure |2l For the Owens model, we set rjp^o = 0.197 Pa s, rjp^oo = 0.003 Pa s, 
rjs = 0.001 Pa s, 07 = 8, and m = 0.75. These parameter values were chosen by Fang 
and Owens 112911 to fit experimental data for the steady-state viscosity of blood, as 
reported by Chien ll^ . The fitted curve and experimental data are reproduced 
in the inset to figure |2l Additionally, Fang and Owens ll^ suggest Ah = 0.004 
s, which leads to Aq = 0.263. In order to compare the two fluid models, we 
assume that the FENE-P model has the same value of 7]p,o and Aq, and that rjs is the 
same. This assumption is based on the expectation that any choice of viscoelastic 
model would have to be compatible with known experimental information on the 
rheology of the fluid, which would, at the least, include a knowledge of the zero 
shear rate viscosity and the relaxation time. Note that for the FENE-P model, while 
there is no necessity to prescribe Aq when the shear rate dependence of viscosity 
is expressed in terms of Wz, it is necessary when represented in terms of y. As is 
well known, the FENE-P model predicts increasing shear thinning with decreasing 
values of the finite extensibility parameter bu- The entire family of curves for the 
FENE-P model shown in figure |2l with values of bu ranging from 100 to 2, does 
not shear thin as rapidly as the Owens model. In particular, it is clear from the 
inset that for the parameters recommended by Fang and Owens 12911 . the FENE-P 
model is unable to capture the rapidity with which blood shear thins, even for 
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Bm = 2. In all cases, in line with expectation, shear thinning first occurs for the 
FENE-P fluid when Wi ~ 0(1). For viscoelastic fluids, the onset of shear thinning 
at Wi ~ 0(1) implies that for shear rates y > 1/Ao, the fluid no longer responds in a 
Newtonian manner, and that the longest time scale for micro-structural relaxation is 
comparable to Aq. Since the Owens model shear thins at Wi « 1, this suggests that 
the characteristic time scale for micro-structural rearrangement is much larger than 
Aq. As mentioned earlier, Aq corresponds to the relaxation time for an aggregate of 
blood cells, which according to Owens |0], are typically of a size that represents the 

— Owens 

greatest proportion of erythrocytes. By defining a Weissenberg number Wi 
for which the Owens model fluid shear thins when it is of 0(1) (see the upper 
horizontal axis of figure IZj, we can estimate that the appropriate relaxation time is 
of order 10“^, which must correspond to much larger structures than a typical blood 
cell aggregate. We do not explore this aspect further here, rather, for the purposes 
of the present paper, we assume that the FENE-P and Owens models are distinct 
constitutive models, which have the same zero shear rate material properties, but 
shear thin significantly differently. As will be discussed in greater detail in the 
sections below, the difference between the models leads to significant differences 
in their behaviour. 

For all the computational results reported here, we set r/p^o, ?]p,oo rjs, di, and m 
at the values recommended by Fang and Owens 12911 . However, we vary Aq (by 
varying Ah) in order to control the inlet Weissenberg number. For the FENE-P 
fluid, we set &m = 100, which is a value commonly used in simulations. As we are 
interested in small blood vessels, we choose the width W of the chaimel to be 100 
pm and Uq = 0.01 m/s, inline with the data reported in Robertson et al. H. The 
value of Re in small blood vessels is well below 1. We have not seen any significant 
difference in the profile shape of the collapsible wall for values of Re in the range 
of 0-1, so we set Re = 0 by setting p = 0. 

Deng and Guidoin 14111 and Intengan et al. 114211 have reported the values of 
Young's modulus (£) for the human artery to be in the range 200-4000 kPa. Zhang 
and co-workers Il43l - l46ll have reported the Young's modulus of the porcine artery 
to be in the range 110-140 kPa, while using two different values for the external 
pressure (8 kPa and 9.3 kPa) in their experimental measurements of the Young's 
modulus. In order to adequately represent the microcirculation, we choose a wide 
range of values for the external pressure pe from 1.2 to 16 kPa and G (which is equal 
to E/3) in the range 30 to 400 kPa. Since Pg = Pe/G, we keep Pg fixed at a constant 
value of 0.04 even though both pg and G are varied. On the other hand, we vary F 
in the range 4.95 X 10“^ to 6.6 x 10“^. For most of the simulations reported here, a 
fixed value of 0.4 W is chosen for the thickness of the solid wall (f), as the artery wall 
thickness to vessel diameter ratio is typically very high in small blood vessels |@]. 
However, in section 13.41 a few results of simulations with varying wall thickness 
are discussed to elucidate the effect of wall thickness. 


3.2. Mesh convergence and the high Weissenberg number problem 


Chakraborty et al. 112211 have established that the flow in a collapsible chaimel 
with a zero-thickness membrane suffers from the high Weissenberg number prob¬ 
lem and have shown that there is a limiting Weissenberg number for each of the 
fluid models beyond which computations fail. Furthermore, this limiting Wi value 
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(a) 




Figure 3: Meshes considered in the current study, (a) Ml, (b) M2, and (c) M3, for t = 0.4W. 


Mesh 

Number of elements 

Number of nodes 

Degrees of freedom 

Ml 

400 

1705 

10972 

M2 

900 

3757 

24072 

M3 

1600 

6609 

42252 


Table 1: Meshes considered in the current study. 


has been shown to increase with mesh refinement. Here we study mesh conver¬ 
gence over a range of parameters for the current geometry, using three different 
meshes M1-M3 for t = 0.41V, as illustrated in figure |3l, with the mesh details given 
in Table [H 

In viscoelastic flow, mesh convergence is generally studied by examining the 
values of the invariants of the conformation dyadic, M. The eigenvalues m; of 
the conformation dyadic represent the square stretch ratios along the principal 
directions of stretching m* for an ensemble of molecules |27, It has been well 
established that the breakdown of viscoelastic computations is typically due to the 
smallest eigenvalue becoming negative in some regions of the flow domain 1221.134 . 


Figure |4] shows the maximum eigenvalue m 3 and minimum eigenvalue mi of 
the conformation tensor as a function of Wi for the Oldroyd-B, FENE-P and Owens 
fluids at F = 4.95 x 10“^ and Pg = 0.04 for t = 0.4W. Figures |4] (a)-(c) clearly exhibit 
the breakdown of viscoelastic computations at a particular value of Weissenberg 
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(c) (f) 


Figure 4: Minimum value of the smallest eigenvalue (mi) and maximum value of the largest 
eigenvalue (m 3 ) in the entire flow domain, for the Oldroyd-B ((a) and (d)), FENE-P ((b) and 
(e)), and Owens model ((c) and (f)), as a function of YJi at E = 4.95 x 10“^ and Pe = 0.04 for 
t = 0.4W. 
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number on each mesh, since the minimum value of mi becomes negative. This 
limiting Wi increases with increase in mesh refinement. 

An increase in Wi leads to a higher maximum m 3 and lower minimum mi 
across the flow domain. While the breakdown value of Wi at each mesh can be 
anticipated from the sudden change of the slope of the curves in the minimum mi 
plots (figure IH (a)-(c)), the curves of the maximum m 3 on various meshes overlap 
with each other (figures |l](d)-(f)). The limiting value of Wi on the M2 mesh for the 
Oldroyd-B, FENE-P and Owens fluids is respectively 0.29,0.38 and 2.13, while the 
mesh converged value of Wi for these models is 0.17, 0.20 and 1.0 respectively. In 
all our analysis, we have ensured that mesh converged values of Wi are used for 
any particular mesh. 

While figure |4] displays the maximum and minimum eigenvalues in the entire 
flow field, figure |5] displays the mean streamwise molecular stretch M^x across 
the channel at the point where the gap between the flexible and rigid walls is a 
minimum, for a fixed value of E, and increasing values of Wi. With increasing 
Weissenberg number, Mxx grows nearly symmetrically from a relatively low value 
in the middle of the gap, to a significantly larger value near the bottom (rigid) 
and top (flexible) walls. Note that in the Oldroyd-B and Owens models, Mxx is 
unbounded due to the infinite extensibility of the Elookean spring in the Elookean 
dumbbell model which underlies these fluid models. Conversely, the existence 
of a upper bound to the mean stretchability of the spring in the FENE-P model 
restricts the maximum value for Mxx, which for bu = 100 is 300. The profiles of Mxx 
for the different fluids in figure |5]clearly reflect this micro-mechanical aspect of the 
models, and confirm that as in other benchmark problems for non-Newtonian flow, 
numerical computations in a 2-D collapsible channel also fail due the development 
of large stresses and stress gradients in certain regions of the flow field, which are 
related to large changes in the conformations of the molecules. 

In their earlier study with a zero-thickness membrane, Chakraborty et al. 11221] 
have shown that the extent of collapse of the membrane also has a significant effect 
on the limiting Weissenberg number. As the gap in the channel becomes narrower 
with decreasing tension in the membrane, the fluid is 'squeezed' leading to a greater 
deformation of the molecules, with a concomitant numerical breakdown at smaller 
values of Wi. One of the parameters that controls the deformability of the wall 
in the current work is T (the other being wall thickness t, whose influence on the 
limiting Wi is discussed subsequently in section |3~^ . Figure [^displays the limiting 
and the mesh converged values of the Weissenberg number for the M2 mesh, as 
a function of E. The limiting Weissenberg number follows the trend Owens > 
FENE-P > Oldroyd-B. For the Oldroyd-B and FENE-P models, the limiting and 
mesh converged Weissenberg numbers appear to increase with an increase in E for 
values of E < 4.0 x 10“^, before levelling off at higher values of E. On the other 
hand, for the Owens model fluid, and Wiconverged increase monotonically 

with E. The Owens model fluid also exhibits the biggest difference between the 
converged and limiting values of Wi. As will be apparent when we discuss the 
shape of the fluid solid interface in figure [71 an increase in T leads to an increase in 
the magnitude of the narrowest channel gap, and consequently an increase in the 
limiting value of Wi. 
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Figure 5: Profile of Mxx across the narrowest channel gap for the Oldroyd-B, FENE-P and 
Owens models, for a range of Weissenberg numbers, at P = 1.98 x 10“^, Pe = 0.04 and 
t = 0.4W. For the FENE-P fluid, we set &m = 100. The distance from the bottom channel is 
scaled by the narrowest gap width Ax/max (see figure 8 (a) for a definition) of the particular 
model. 
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Figure 6: Maximum mesh converged value of Wi and the limiting Wi, for the three fluid models, 
for computations carried out with the M2 mesh, at Pe = 0.04 and t = 0.4 W, as a function of F. 


3.3. Interface shape and velocity contours 

Figure IZl explores the deformation of the finite-thickness solid wall, while inter¬ 
acting with the different fluids, at a fixed value of t = 0.4W. While figures [71 (a)- 
(c) investigate the shape of the fluid-solid interface for different values of F at 
Wi = 0.1, figures [Zl(d)-(f) examine the dependence of the interface profile on Wz for 
F = 4.95 X lO""^. The extraordinary variation in the shape of the elastic solid with 
varying elasticity parameter F is immediately apparent from figures [71 (a)-(c). In 
particular, there is a stark contrast in the response of the solid to the flow of differ¬ 
ent viscoelastic fluids. Except in the case of the Owens model, increasing F leads 
to a movement of the deformable solid from being within the charmel (concave 
downwards) to bulging out of the charmel (convex upwards) due to action of the 
forces exerted by the flowing fluid. At the relatively low value of Wz = 0.1 there is 
no discernible difference between the Newtonian, Oldroyd-B and FENE-P fluids. 
For the Owens model on the other hand, while the elastic solid remains concave 
downwards for the entire range of F values, there is a decrease in the minimum 
channel gap with increasing F. This behaviour is related to the significant differ¬ 
ence in the force field generated in the Owens model fluid due to flow, as discussed 
shortly. 

It is appropriate to note here that in our earlier investigation of viscoelastic flow 
in a 2D charmel with a zero-thickness membrane [[220, the fluid-solid interface was 
always observed to be concave downwards for all the viscoelastic fluids, at all val¬ 
ues of membrane tension. Indeed, in contrast to the situation for a finite thickness 
solid, with decreasing tension, the zero-thickness membrane moves further into 
the charmel, with a concomitant decrease in the narrowest charmel gap. 
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Figure 7: The shape of the fluid-solid interface in a 2-D collapsible channel for the Oldroyd- 
B ((a) and (d)), FENE-P ((b) and (e)) and Owens models ((c) and (f)), compared with the 
profile for a Newtonian fluid. Note that Wi is 0.1 in (a)-(c) and E is 4.95 x 10“^ in (d)-(f). In 
(a)-(c) different symbols represent different values of T (□: 4.95 x 10“^, O; 1.98 x 10“^, 'A' : 
3.0 X 10“^, +: 3.96 x 10“^, X: 4.95 x 10“^ and A: 6.6 x 10“^). Lines with the same colour as 
the symbols represent the predictions of a Newtonian fluid for identical values of E. 
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Figure 8: (a) Schematic diagram defining the position of maximum deformation 

{Axjnax, Syntax) ■ (b) Dependence of {Axmax, Aymax) on r at a fixed value of Wi = 0.1, and 
(c) on Wi at a fixed value of F = 4.95 x 10“"^. In (b) and (c) the arrows indicate the direction 
of increasing F and Wi, respectively. The range of Wz for the Oldroyd-B, FENE-P and 
Owens models are 0.01-1.508, 0.01-2.372 and 0.01-7.9, respectively. 
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Figure 9: Contours of axial velocity (Vx) in the flow domain, for Newtonian (black), Oldroyd- 
B (red), FENE-P (blue) and Owens (green) fluids at Pg = 0.04, t = 0.4W, for two different 
values of Weissenberg number Wi = 0.1 ((a) and (c)) and Wi = 0.5 ((b) and (d)). Note that 
r = 1.98 X 10“^ in (a)-(b) and E = 4.95 x 10“^ in (c)-(d). The upper boundary in these figures 
reflects the shape of the interface at the corresponding parameter values. 


At a fixed value of elasticity parameter T, while figure[Zl(f) indicates that Wi has 
no effect on the shape of the deformable solid in the case of the Owens model (which 
remains concave downwards), it has a noticeably different effect for the Oldroyd-B 
and FENE-P fluids. Both fluids cause the elastic solid to bulge outwards. However, 
the extent of this bulge decreases more rapidly for the FENE-P fluid with increasing 
Wi. In the case of the Owens model, at these values of Wi, shear thiiming is nearly 
complete, and there is consequently no change discernible in the membrane shape. 
On the other hand, the onset of shear thinning for the FENE-P model is responsible 
for the observed variation in the predicted membrane shape from that for an 
Oldroyd-B fluid. 

Figure [71 indicates that the deformation of the solid wall occurs along both 
axial and vertical directions for all the fluid models, with the extent of movement 
depending on the values of E and Wi. By defining the position of maximum 
deformation as the point on the elastic solid furthest in the vertical direction from 
the horizontal surface, this dependence can be examined more systematically. The 
precise location of the position of maximum deformation is given by the co-ordinate 
pair {Axmax, Ay^a^, as shown schematically in figure |8] (a), which measures the 
maximum displacement from the centre of the elastic solid when it is horizontal. 
Figures|8](b) and|8](c) track the position of maximum deformation for varying E and 
Wi, and correspond to the set of figures [Zl(a) to (c) and[Zl(d) to (f), respectively. The 
movement of the elastic solid from being concave downwards to convex upwards 
in the case of varying E, and the downward movement with increasing Wi are 
clearly captured in this description. The relative immobility of the solid in the case 
of a flowing Owens model fluid is also clearly revealed. 

Figure |9| compares the velocity contours predicted by a Newtonian fluid with 
those of an Oldroyd-B, FENE-P and Owens' fluid at Pg = 0.04 and t = 0.4W, for 
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two different values of T and Wi. The qualitative and nearly quantitative similarity 
of the velocity contours for the Oldroyd-B and FENE-P models with those for a 
Newtonian fluid at Wi = 0.1 and Wi = 0.5 suggests that fluid rheology does not 
have a significant influence on the velocity field at these values of Wi. On the other 
hand, the dominant influence appears to be the value of the elasticity parameter E, 
which determines the shape of the fluid-solid interface. While the velocity contours 
for all the fluids are qualitatively similar to each other when the elastic solid lies 
within the charmel (figures |9] (a) and (b) for the Oldroyd-B and FENE-P models, 
and figures |9] (a) to (d) for the Owens model), there is a qualitative change when 
the elastic solid lies outside the channel. The latter situation occurs only for the 
Oldroyd-B and FENE-P models for E = 4.95 x 10“^, as displayed in figures |9](c) and 
(d). 

3.4. Influence of wall thickness on interface shape and limiting Weissenberg number 

The wall thickness t does not appear directly in the governing equations, or 
in the boundary conditions, rather it determines the size of the solid domain. 
We can anticipate, however, that its influence on the interface shape and limiting 
Weissenberg number will be similar to that of the elasticity parameter T. In very 
simple terms, increasing T for a given state of stress at the fluid-solid interface due 
to the flow of viscoelatic fluid in the charmel, leads to a larger strain in the solid, 
since (with all other parameter values fixed) an increase in T implies a decrease 
in the shear modulus of the solid, G. Similarly, for a given state of stress at the 
fluid-solid interface, decreasing the wall thickness t would lead to a larger strain in 
the solid, since there is "less" solid material over which to distribute the resultant 
stress in the solid. 

These arguments are home out in figures [TOl, which explores the deformation 
of the finite-thickness solid wall, while interacting with the different fluids, at a 
fixed value of T = 4.95 x 10“"^. While figures [TOl (a)-(c) investigate the shape of the 
fluid-solid interface for different values of t at Wi = 0.1, figures [TOl (d)-(f) examine 
the dependence of the interface profile on Wi for f = O.IW. 

Figures [To] (a) and (b) indicate that for E = 4.95 x 10““^, at all values of t examined 
here, the solid wall bulges outward, and there is no discernible difference between 
the Oldroyd-B and FENE-P models, and a Newtonian fluid. For all these fluids, 
increasing t leads to a decrease in the extent of deformation of the solid wall. 
The profiles are similar to those observed previously in figures [71 (a) and (b) for 
decreasing values of E. In the case of the Owens model fluid, the fluid-solid 
interface lies within the charmel for all the examined values of wall thickness. 
However, there is a distinct change observed in the shape of the interface, with the 
profile becoming more symmetric as the value of t increases. 

The change in interface shape with increasing Wi, at fixed values of t and 
T, displayed in figures [T0[ (d)-(f) is very similar to that observed previously in 
figures [T] (d)-(f), with the extent of shear thinning playing the dominant role in 
determining the shape. The essential difference between the two sets of figures 
appears to be in the loss of symmetry in the interface profile for the smaller value 
of t. 

As discussed earlier, both the mesh converged and limiting Weissenberg num¬ 
bers for the Owens model increase with increasing T because of an increase in the 
narrowest charmel gap (see figures [6[ and figure 0(c)). This argument is consistent 
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Figure 10: The shape of the fluid-solid interface for the Oldroyd-B ((a) and (d)), FENE-P ((b) 
and (e)) and the Owens model fluids ((c) and (f)), compared with the profile for a Newtonian 
fluid, at different values of wall thickness t and Wi, at a fixed value of E = 4.95 x 10“'^. The 
value of Wi is fixed at 0.1 in (a)-(c), while t is fixed at O.IW in (d)-(f). In (a)-(c) different 
symbols represent different values of t (★ : O.IW, O; 0.2W, □: 0.4W). Lines with the same 
colour as the symbols represent the predictions of a Newtonian fluid for identical values 
of t. 
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Figure 11: Limiting Weissenberg number for the three fluid models, for computations carried out 
with the M2 mesh, at Pe — 0.04 and F = 4.95 x 10“^, as a function of the non-dimensional wall 
thickness f/W. 


with the behaviour of Wiumiting displayed in figure [H] for the Owens model, for 
increasing values of t/W (see also figure [lO] (c)). In both these situations corre¬ 
sponding to the Owens model, the elastic wall always lies within the chaimel for 
all the values of the various parameters considered here. We have seen previously 
in figurelH for both the Oldroyd-B and FENE-P models, that Wiumiting and Wiconverged 
are not sensitive to changes in E for values of E that correspond to the situation 
where the elastic wall lies outside the charmel. On the other hand, figure [H] sug¬ 
gests that even though the elastic wall lies outside the charmel for E = 4.95 x 10“^ 
(see figures HOja) and (b)), Wiumiting increases with t/W, until t/W = 0.2, before 
levelling off. 


3.5. Pressure and stresses 


Patankar et al. Il49l] have shown analytically that for any constitutive model of 
the form. 


fljD -F a2D -i- a;^P -i- aPT — 0 


( 20 ) 


where, fli, a 2 , ^ 3 , and ^4 are scalar functions of the invariants of D and T, and D 

V 

and T are the upper convected time derivatives of D and T, the normal component 
of extra stress on a rigid body surface will be zero. Chakraborty et al. [1221] have 
shown numerically that this is true even in the case of flow past a deformable 
zero-thickness membrane, for all the three viscoelastic fluids considered here. In 
the present instance as well, we find that the normal component of stress on the 
elastic wall is solely due to pressure. 

Figure [12] examines the effect of E and Wi on the non-dimensional pressure P 
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Figure 12: Dependence of the pressure profile along the flexible membrane on Wi and F, 
for the Oldroyd-B ((a) and (d)), FENE-P ((b) and (e)) and Owens models ((c) and (f)), 
respectively. The lines in (a)-(c) are for a Newtonian fluid. Note that E = 1.98 x 10“"^ in 
(d)-(f) and Wi = 0.1 in (a)-(c). 
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exerted by the different fluids on the elastic solid. At Wz = 0.1, the increase in 
P with increasing T for the Oldroyd-B and FENE-P models is nearly identical to 
that for a Newtonian fluid, as can be seen from figures [l2l (a) and (b). Notably, 
for all these fluids, a distinct change occurs in the shape of the pressure profile 
for E > 1.98 X 10“"^. For values of E less than or equal to this value, the pressure 
profile is relatively constant upstream of the position of maximum deformation, 
before decreasing relatively rapidly downstream to a constant value. (The vertical 
lines in the figure denote the x-position of maximum deformation, with the colour 
coordinated to match the corresponding E value). On the other hand, for values 
of E > 1.98 X 10“^, the decrease in pressure from the location where the fluid 
flows under the deformable solid to the location where it exits, is much more 
uniform. As can be seen from figures [71 (a) and (b), the change in the shape of the 
pressure profile is correlated with the change in interface shape that occurs around 
E ~ 3 X 10“"^, which is approximately the value at which the elastic solid moves 
from being concave downwards within the chaimel to bulging outwards from the 
channel. In the case of the Owens model, even though the pressure increases with 
increasing E, the shape of the pressure profile remains unchanged, since the elastic 
solid is always concave downwards in shape (see figures 1121 (cl and [71 (c)). Another 
notable aspect is that the magnitude of pressure at any point along the interface is 
significantly lower for the Owens model compared to that for all the other fluids. 
This can be attributed to the significant decrease in viscosity that occurs for the 
Owens model fluid when it flows under the deformable solid. 

Figure [12] (d)-(e) display the effect of Wz on the pressure profile for a fixed value 
of E = 4.95 X 10“"^. At this value of E, as seen earlier in figures [7[ (d)-(f), for all the 
values of Wz considered here, the elastic solid bulges outwards from the chaimel 
due to interaction with the Oldroyd-B and FENE-P fluids, while it is concave 
downwards for the flow of an Owens model fluid. In the former two cases, with 
increasing Wz, there is a clear decrease in the pressure that the fluid exerts on the 
downstream end of the elastic solid, with the decrease being more substantial for 
the FENE-P fluid. This correlates with the decrease in the bulge of the elastic solid 
seen earlier in figures 0 (d) and (e). For the Oldroyd-B fluid, there also appears to 
be a slight increase in pressure at the upstream end of the elastic solid. In the case 
of the Owens model fluid, neither the interface shape nor the pressure profile are 
significantly altered by the variation in Wz. 

Figures [T3[ display the pressure contours under the collapsible wall for all the 
fluids at Pe = 0.04 and t = 0.4W, for different values of E and Wz. As was ob¬ 
served earlier in the case of velocity contours in figures [U the qualitative shape 
of the pressure contour depends strongly on whether the elastic solid lies within 
or outside the channel. Thus, shapes for all the fluids are similar to each other in 
figures [T3[ (a) and (b), while the shapes for the Newtonian, Oldroyd-B and FENE-P 
fluids differ qualitatively from that of Owens model fluid in figures [T3[ (c) and (d), 
since the elastic wall lies inside the charmel in the latter case while lying outside 
for the former. 

The pressure contours for the Oldroyd-B and FENE-P fluids are quantitatively 
similar to those for the Newtonian fluid only for the case when Wz = 0.1 and 
E = 1.98 X 10“"^ (figure [T3[ (ah. On the other hand, they depart from the Newtonian 
contours at the remaining values of Wz and E displayed in figures [13] (b) to (d). 
This is consistent with the behaviour of the pressure profiles along the fluid-solid 


22 



Figure 13: Contours of pressure in the flow domain, for Newtonian (black), Oldroyd-B 
(red), FENE-P (blue) and Owens (green) fluids at Pe = 0.04, t = 0.4W for two different 
values of Weissenberg number Wi = 0.1 ((a) and (c)) and Wi = 0.5 ((b) and (d)). Note that 
E = 1.98 X 10“^ in (a)-(b) and E = 4.95 x 10“^ in (c)-(d). 


interface observed in figures IT^ (a), (b), (d) and (e). 

A different perspective on fluid pressure in the channel is provided in figure [HI 
where the pressure drop AP in the channel between the entrance and exit to the 
region beneath the elastic solid, is displayed. As seen earlier in figures |Z1 (a)-(b), 
with increasing E, the interface shape for the Oldroyd-B and FENE-P fluids moves 
from being concave downwards to convex upwards. Figure [14] (a) shows that this 
is accompanied by an increase in AP. Interestingly, the rate of change of AP with 
E has a point of inflection around E ~ 3 X 10“^, which is approximately the value 
at which the elastic solid becomes horizontal (see figures |ZI (a)-(b), and inset to 
figure IT^ (a)). 

A striking manifestation of differences in the prediction of a macroscopic prop¬ 
erty, because of differences in fluid rheology, is displayed in figure [HI (b), where 
the dependence of pressure drop AP on Weissenberg number Wi is plotted. The 
Owens model fluid has a nearly constant pressure drop because the fluid has un¬ 
dergone significant shear thinning, and has an almost constant viscosity under the 
deformable elastic solid for all values of Wi. For the Oldroyd-B model on the other 
hand, which is a constant viscosity fluid, there appears to be a very slight increase 
in AP. Clearly, the decrease in pressure at the downstream end of the chaimel, is 
more than made up with the increase at the upstream end. For the FENE-P fluid, 
the increasing shear thinning with increasing Wi is reflected in figure [H (b) with 
the observed decrease in AP. 

Figure [15] explores the dependence of the axial component of the conformation 
tensor M^x, along the flexible wall, on E and Wi. Figures [15[ (a)-(c) show that an 
increase in E leads to a decrease in the degree of stretching experienced by the 
micro-structural elements, and that the value of Mxx in the Owens model is much 
less sensitive to the value of E compared to the Oldroyd-B and FENE-P fluids. 
For the latter two fluids, for values of E < 3 x 10“^, the elastic solid is concave 
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Figure 14: Dependence of pressure drop AP in the channel for the Oldroyd-B, FENE-P and 
Owens models on (a) E at a fixed value of Wz=0.1 and (b) Wi at E = 4.95 x 10“^. Note that 
for a Newtonian fluid, AP = 0.1 in (b)micro-structural. The curves terminate at the limiting 
Weissenberg number for each model. 
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Figure 15: Dependence of the axial component of the conformation tensor Mxx on Y, for (a) 
Oldroyd-B, (b) FENE-P, and (c) Owens models, at Wi = 0.1, and dependence of Mxj on Wi, 
for (d) Oldroyd-B, (e) EENE-P, and (f) Owens models, at E = 4.95 x 10“^. 
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downwards. As a result, the Mxx profile has a maximum at the location in the 
channel where the gap is narrowest. As the elastic solid moves out of the channel, 
there is a significant relaxation in the degree to which the micro-structural elements 
are stretched. 

The correlation between interface shape and Mxx profile is more strikingly re¬ 
vealed in figures [15] (d)-(e), where the dependence of Mxx on Wi is explored at a 
constant value T = 4.95 x 10“^. Since the interface always bulges outwards for the 
Oldroyd-B and FENE-P fluids at this value of E, the highest stretch occurs at the 
inlet and outlet to the deformable region, in contrast to the situation for the Owens 
model, where the elastic solid is always concave downwards, and consequently, 
the maximum stretch is always at the location of the narrowest gap. At high values 
of IVi, the shear thinning experienced by the FENE-P fluid appears to lead to a 
more uniform stretching along the length of the channel. In all cases, however, as 
might be anticipated, an increase in Wz leads to an increase in stretching. 

Figures [16] (a)-(c) show the contour plots of the mean streamwise molecular 
stretch Mxx for the Oldroyd-B, FENE-P and Owens models, at E = 1.98 x lO""^ and 
Wz = 0.1. As can be seen from figures [15] (a)-(c), these parameter values correspond 
to the situation where the elastic wall lies within the channel for all the fluids, 
and the values of Mxx for the Oldroyd-B and FENE-P models are close to each 
other along the entire length of the channel. Further, for all the models, the largest 
value of Mxx occurs below the collapsible wall at the minimum gap location. All 
these observations are clearly reflected in figures [16] (a)-(c), both in the shape of the 
interface and in the values of Mxx corresponding to the various contour lines. 

Finally, the dependence of the total shear stress on the elastic solid, r® -i- on 
the parameters E and Wz, is examined in figure [17] for the three viscoelastic fluids. 
Once again, there is close parallel between the shape of the fluid-solid interface 
and the shear stress on the wall. Indeed, the shear stress profiles are either concave 
downwards or convex upwards in complete synchrony with the interface shape. 
In contrast to the zero-thickness membrane model, where the shear stress on the 
membrane has no influence on membrane shape because of the use of a boundary 
condition that only accounts for the influence of the normal stress, in the present 
model, both the pressure and the shear stress are responsible for the membrane 
shape. As a result, a much greater variety of interface shapes is observed for a 
finite thickness elastic solid. 

4. Conclusions 

We have introduced a new geometry, whose central feature is the existence of 
fluid-structure interaction, into the lexicon of standard benchmark non-Newtonian 
flow computations. The role that the presence of a deformable membrane plays 
in the development of a complex flow field in the channel has been examined, 
and the relationship of the upper limit to the Weissenberg number to molecular 
conformations at various locations in the flow domain, has been delineated. The 
shape of the membrane as a function of a membrane elasticity parameter E, and of 
the Weissenberg number Wz has been studied, and the change in shape has been 
used as an indication of the extent of fluid-structure interaction. The nature of the 
coupling between macroscopic observables such as velocity, stress and conforma¬ 
tion fields, and various rheological features of the three viscoelastic fluid models 
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Figure 16: Contours of Mxx in the flow domain for the Oldroyd-B, FENE-P and Owens 
models at E = 1.98 x 10“^, Pe = 0.04, t = 0.4W and Wi = 0.1. 


27 






















' 

° □ 

1 

1 

□ 

□ 

1 

□ 

^ □ 

1 

1 

° O 

1 

;★ 

'A 

1 

o 

o 

□ 

A 

□ 

--A-< 

A 

' 

A 

A o 

(I 

1 

*■ . 
o ^ 

★ 

□ r=4.95xl0'^ 

6 

1 

1 

1 

1 ★ 
0| ^ 

1 

o r=1.98xl0“'‘ 

1 

1 

1 

1 



1 Oldroyd-B 

^ r=3.96xl0 

1 

1 

★ r=4.95xl0“'^ 

1 

1 

1 

1 


.5 -1.25 0 1.25 2.5 


xAV 

(a) 


° □ 

1 1 1 

1 1 1 

□ 

□ 

! ! ! □ 


^ □ 

1 1 1 


o ★ 

O ★ a, 

4 ^ 

' ' ' □ 

; V 

\ JA O 

nil A . 

o 

A O 

A. 

: :: o 

A 

★ 

□ r=4.95xl0'^ 
o r=1.98xl0“'‘ 

O ' ' 

' Pi ° 

1 1 1 

1 1 1 

1 1 1 

★ 


1 1 1 FENE-P 


^ r=3.96x10 
★ r=4.95xl0“'‘ 

1 1 1 

1 1 1 

1 1 1 

1 1 1 



,.5 -1.25 0 1.25 2.5 


xAV 

(b) 


A ®aB o T 

= '8 8*^ 

1 

i;“ A 

O □' 

: o ★ 

★ 

o 

o 

" A 

4 

1 ★ 

□ r=4.95xl0'^ 


O r=l.98x10“'^ i 

" Owens 

^ r=3.96xl0“'‘ 


★ r=4.95xl0“'‘ : 



:.5 -1.25 0 1.25 2.5 

xAV 


(c) 



(d) 



(e) 



(f) 


Figure 17: Dependence of the tangential component of stress on F, for (a) Oldroyd-B, 

(b) FENE-P, and (c) Owens models, at Wi = 0.1, and dependence of r® + tJ’ on Wi, for (d) 
Oldroyd-B, (e) FENE-P, and (f) Owens models, at E = 4.95 x 10“^. 
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used in this study, has been explored in some detail. 

There are many aspects of viscoelastic flows in two dimensional collapsible 
channels that remain to be studied, (i) The use of a constitutive model that accounts 
for thixotropy is an important feature, since the aggregation of blood cells in regions 
of low shear rate can lead to rheological properties that depend locally on micro- 
structural dynamics. Owens model in its most general form does account for 
thixotropy 101, and as mentioned earlier, lolov et al. ll^ have recently developed 
a finite element method for solving the Owens model in its complete generality, 
(ii) Even though there exists an upper limit to the Weissenberg number at which 
computations fail for each mesh, we have not encountered, in our admittedly 
limited simulations, a situation where this upper limit has not changed in spite of 
mesh refinement. It would be interesting to see if the use of a log-conformation 
tensor formalism leads to much higher upper limits to the Weissenberg number 
for all the models, (iii) The multiple modes of instabilities that arise for flow 
in collapsible chaimels, and the rich behaviour that occurs in unsteady flows, has 
been extensively investigated for Newtonian fluids. We hope that the present work 
provides a starting point for similar studies in the context of viscoelastic fluids. 
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